function [TotalTime, mvcount, TotalBytes, L2abserrpower, L2relerrpower] = ...
  testnorm(N, nLevel)
if(1)

  h = (1/N);
  V = 1 + rand(N,N);
  % V = ones(N,N);

  global Tree;
  disp('Factorization');
  Tree = five_setup(V,N,h);
  disp('Extraction');
  DiagExact = five_extract(N,Tree);
end


if(1)
  svdfact = 1;
  epsilon = 1e-6;


  global mvcount TotalTime;
  mvcount = 0;
  TotalTime = 0;
  [DiagVec, DiagBlock, nDiagBlock, SampleList, RKTree] = ...
    ConstructHMatrix( N, nLevel, epsilon, svdfact );
  disp(['Number of matrix-vector applications is ' num2str(mvcount)])
  fprintf('\n')

  disp(['relative L^1 error of the diagonal elements']);
  diagrelerrL1 = sum(sum(abs(DiagVec - DiagExact))) ./ ...
    sum(sum(abs(DiagExact)))

  disp(['absolute Max error of the diagonal elements']);
  diagabserrMax = max(max(abs(DiagVec - DiagExact)))

  disp(['relative Max error of the diagonal elements']);
  diagrelerrMax = max(max(abs(DiagVec - DiagExact))) ./ ...
    max(max(abs(DiagExact)))

end

% Estimate 2-norm using random method
if(0)
  nSample = 100;
  u = rand(N,N,nSample);
  v1 = WholeTreeApply( ...
    RKTree, DiagBlock, u, N, nSample, nLevel+1, nDiagBlock);
  v2 = KernelApply( N, reshape(u, N*N, nSample) );
  disp(['Estimated L^2 error using ' num2str(nSample) ' tests']);
  L2err = max( sqrt(sum((v1-v2).^2,1)) ./ sqrt(sum(v1.^2,1)) )
end

% Estimate 2-norm using power method
if(0)
  nSample = 1;
  nPowerMax = 500;
  u = rand(N,N,nSample);
  v = reshape(u, N*N, nSample);
  for iter = 1 : nPowerMax
    if( mod(iter, 20) == 0 )
      fprintf('%5i th power iteration for (G-G_approx)\n\n', iter);
    end
    w = reshape(v, N, N, nSample);
    v1 = WholeTreeApply( ...
      RKTree, DiagBlock, w, N, nSample, nLevel+1, nDiagBlock);
    v2 = KernelApply( N, reshape(w, N*N, nSample) );
    vdiff = reshape(v1-v2, N*N, nSample);
    vt = vdiff(:,1)./norm(vdiff(:,1))./v(:,1);
    if( std(vt) < 1e-2 )
      fprintf('%5i iterations to reach std < 1e-2 \n\n', iter);
      break;
    end
    v = vdiff;
    for j = 1 : nSample
      v(:,j) = v(:,j) / norm(v(:,j),2);
    end 
  end
  w = reshape(v, N, N, nSample);
  v1 = WholeTreeApply( ...
    RKTree, DiagBlock, w, N, nSample, nLevel+1, nDiagBlock);
  v2 = KernelApply( N, reshape(w, N*N, nSample) );
  vdiff = reshape(v1-v2, N*N, nSample);
  L2abserrpower = abs(sum(v.* vdiff, 1));
  disp(['Estimated absolute L^2 error using power method']);
  L2abserrpower
end

% Estimate the L^2 norm of G matrix to calculate relative L^2 error
if(0)
  nSample = 1;
  nPowerMax = 500;
  u = rand(N,N,nSample);
  v = reshape(u, N*N, nSample);
  for iter = 1 : nPowerMax
    if( mod(iter, 20) == 0 )
      fprintf('%5i th power iteration for (G-G_approx)\n\n', iter);
    end
    w = v;
    v1 = KernelApply( N, reshape(w, N*N, nSample) );
    vt = v1(:,1)./norm(v1(:,1))./w(:,1);
    if( std(vt) < 1e-2 )
      fprintf('%5i iterations to reach std < 1e-2 \n\n', iter);
      break;
    end
    v = v1;
    for j = 1 : nSample
      v(:,j) = v(:,j) / norm(v(:,j),2);
    end 
  end
  v1 = KernelApply( N, reshape(v, N*N, nSample) );
  L2G = sum(v1.* v, 1);
  disp(['Estimated relative L^2 error using power method']);
  L2relerrpower = L2abserrpower / L2G
end



% Estimate 2-norm using power method. The criterion is to measure the
% eigenvalue difference.
if(1)
  nSample = 1;
  nPowerMax = 500;
  u = rand(N,N,nSample);
  v = reshape(u, N*N, nSample);
  Eold = 1;
  Enew = 0;
  for iter = 1 : nPowerMax
    if( mod(iter, 20) == 0 )
      fprintf('%5i th power iteration for (G-G_approx)\n\n', iter);
    end
    w = reshape(v, N, N, nSample);
    v1 = WholeTreeApply( ...
      RKTree, DiagBlock, w, N, nSample, nLevel+1, nDiagBlock);
    v2 = KernelApply( N, reshape(w, N*N, nSample) );
    vdiff = reshape(v1-v2, N*N, nSample);
    Enew = abs(sum(v.* vdiff, 1));
    if( abs(Enew - Eold)./ Enew < 1e-3 )
      fprintf('%5i iterations to reach abs(Enew - Eold)./ Enew < 1e-3 \n\n', iter);
      break;
    end
    v = vdiff;
    Eold = Enew;
    for j = 1 : nSample
      v(:,j) = v(:,j) / norm(v(:,j),2);
    end 
  end
  L2abserrpower = Enew
  disp(['Estimated absolute L^2 error using power method']);
  L2abserrpower
end

% Estimate the L^2 norm of G matrix to calculate relative L^2 error
if(1)
  nSample = 1;
  nPowerMax = 500;
  u = rand(N,N,nSample);
  v = reshape(u, N*N, nSample);
  Eold = 1;
  Enew = 0;
  for iter = 1 : nPowerMax
    if( mod(iter, 20) == 0 )
      fprintf('%5i th power iteration for (G-G_approx)\n\n', iter);
    end
    w = v;
    v1 = KernelApply( N, reshape(w, N*N, nSample) );
    Enew = abs(sum(v1.*v,1));
    if( abs(Enew - Eold)./ Enew < 1e-3 )
      fprintf('%5i iterations to reach abs(Enew - Eold)./ Enew < 1e-3 \n\n', iter);
      break;
    end
    Eold = Enew;
    v = v1;
    for j = 1 : nSample
      v(:,j) = v(:,j) / norm(v(:,j),2);
    end 
  end
  L2G = Enew;
  disp(['Estimated relative L^2 error using power method']);
  L2relerrpower = L2abserrpower / L2G
end




% Estimate memory cost for both the low rank tree and the diagonal and
% near-off-diagonal blocks. We only count the real numbers. The indices
% are not counted.
%
if(1)
  TotalBytes = 0;
  for iLevel = 1 : nLevel
    S = RKTree{iLevel}.RKBlockU;
    infoS = whos('S');
    TotalBytes = TotalBytes + infoS.bytes;
    S = RKTree{iLevel}.RKBlockV;
    infoS = whos('S');
    TotalBytes = TotalBytes + infoS.bytes;
  end
  for iDiagBlock = 1 : nDiagBlock
    S = DiagBlock{iDiagBlock}.Mat;
    infoS = whos('S');
    TotalBytes = TotalBytes + infoS.bytes;
  end
  fprintf('Total Memory Cost for RKTree and DiagBlock (real numbers only) is %6.3f MB\n', ...
    TotalBytes/1024^2);
end
